d <-#read_csv(paste0(home, "/data/for-analysis/dat.csv")) |> read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv") |>filter(age_ring =="Y", # use only length-at-age by filtering on age_ring!area %in%c("VN", "TH"))
Rows: 338460 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): age_ring, area, gear, ID, sex, keep
dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 43,886 rows (13%), 294,574 rows remaining
# Sample size by area and cohortns <- d |>group_by(cohort, area) |>summarise(n =n())
group_by: 2 grouping variables (cohort, area)
summarise: now 431 rows and 3 columns, one group variable remaining (cohort)
# Minimum number of observations per area and cohortd$area_cohort <-as.factor(paste(d$area, d$cohort))d <- d |>group_by(area_cohort) |>filter(n() >100)
# Minimum number of observations per area, cohort, and aged$area_cohort_age <-as.factor(paste(d$area, d$cohort, d$age_bc))d <- d |>group_by(area_cohort_age) |>filter(n() >10)
# We can also consider removing individuals where the SE of k is larger than the fitIVBG |>#mutate(keep = ifelse(k > quantile(k_se, probs = 0.95), "N", "Y")) |>mutate(keep =ifelse(k < k_se, "N", "Y")) |>#filter(row_id < 10000) |>ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr, color = keep)) +geom_point(alpha =0.5) +facet_wrap(~keep, ncol =1) +geom_errorbar(alpha =0.5) +NULL
mutate: new variable 'keep' (character) with 2 unique values and 0% NA
# Add back cohort and area variablesIVBG <- IVBG |>left_join(d |>select(ID, area_cohort)) |>separate(area_cohort, into =c("area", "cohort"), remove =TRUE, sep =" ") |>mutate(cohort =as.numeric(cohort))
Joining with `by = join_by(ID)`
left_join: added 2 columns (area_cohort_age, area_cohort)
> rows only in x 0
> rows only in y (146,393)
> matched rows 142,023 (includes duplicates)
> =========
> rows total 142,023
mutate: converted 'cohort' from character to double (0 new NA)
# Summarise and save for sample sizesample_size <- IVBG |>group_by(area) |>summarise(n_cohort =length(unique(cohort)),min_cohort =min(cohort),max_cohort =min(cohort),n_individuals =length(unique(ID)),n_data_points =n())
group_by: one grouping variable (area)
summarise: now 10 rows and 6 columns, ungrouped
write_csv(sample_size, paste0(home, "/output/sample_size.csv"))# Compare how the means and quantiles differ depending on this filteringIVBG_filter <- IVBG |>drop_na(k_se) |>#filter(k_se < quantile(k_se, probs = 0.95)) |> filter(k_se < k)
group_by: 2 grouping variables (cohort, area)
summarize: now 306 rows and 6 columns, one group variable remaining (cohort)
ungroup: no grouping variables
# No difference at all really. We should probably just discuss that with this model, achieving biologically reasonable parameter values and a good fit to data are sometimes two different things. In our case, we just want a representative value of the growth (as in time to reach average maximum size in the population) that accounts for the entire growth trajectory of an individual, for each area and cohort.
BT 1996 has a very high k, so I’ll inspect it in more detail
Rows: 380 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): area, model
dbl (2): year, temp
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (cohort)
VBG <- VBG |>left_join(pred_temp, by =c("area", "cohort"))
left_join: added 2 columns (temp, model)
> rows only in x 0
> rows only in y ( 74)
> matched rows 306
> =====
> rows total 306
# Save data for map-plotcohort_sample_size <- IVBG |>group_by(area, cohort) |>summarise(n =n()) # individuals, not samples!
group_by: 2 grouping variables (area, cohort)
summarise: now 306 rows and 3 columns, one group variable remaining (area)
VBG <-left_join(VBG, cohort_sample_size, by =c("cohort", "area"))
left_join: added one column (n)
> rows only in x 0
> rows only in y ( 0)
> matched rows 306
> =====
> rows total 306
write_csv(VBG, paste0(home, "/output/vbg.csv"))# Calculate the plotting order (also used for map plot)order <- VBG |>ungroup() |>group_by(area) |>summarise(min_temp =min(temp)) |>arrange(desc(min_temp))
ungroup: no grouping variables
group_by: one grouping variable (area)
summarise: now 10 rows and 2 columns, ungrouped
order
# A tibble: 10 × 2
area min_temp
<chr> <dbl>
1 SI_HA 7.83
2 BS 6.22
3 BT 5.87
4 FM 5.87
5 SI_EK 5.49
6 MU 5.16
7 FB 5.04
8 HO 3.99
9 JM 3.37
10 RA 3.06
write_csv(order, paste0(home, "/output/ranked_temps.csv"))nareas <-length(unique(order$area)) +2# to skip the brightest colors that are hard to readcolors <-colorRampPalette(brewer.pal(name ="RdYlBu", n =10))(nareas)[-c(6,7)]
Plot VBGE fits
# Sample 30 IDs per area and plot their data and VBGE fitsset.seed(4)ids <- IVBG |>distinct(ID, .keep_all =TRUE) |>group_by(area) |>slice_sample(n =30)
ggsave(paste0(home, "/figures/vb_pars.pdf" ), width =17, height =17, units ="cm")
Fit Sharpe-Schoolfield model to k
By area
model <-'sharpeschoolhigh_1981'# Get starting values on full dataset for Sharpe-Schoolfield modeldat <- VBG |>select(k_median, temp, area) |>rename(rate = k_median)
pivot_wider: reorganized (parameter, Estimate) into (r_tref, e, eh, th) [was 4x2, now 1x4]
# Make predictions on new datanew_data_all <-data.frame(temp =seq(min(dat$temp), max(dat$temp), length.out =100))pred_all <-augment(fit_all, newdata = new_data_all) |>mutate(area ="all")
mutate: new variable 'area' (character) with one unique value and 0% NA
# To get confidence intervals around our nls fits, we can do use "investr"# https://stackoverflow.com/questions/61341287/how-to-calculate-confidence-intervals-for-nonlinear-least-squares-in-r# https://stackoverflow.com/questions/52973626/how-to-calculate-95-prediction-interval-from-nlspred_all$upr <-as_tibble(predFit(fit_all, newdata = new_data_all, interval ="confidence", level =0.95))$uprpred_all$lwr <-as_tibble(predFit(fit_all, newdata = new_data_all, interval ="confidence", level =0.95))$lwr# Add t_opt (not correct equation I think!) from Padfield 2021 ISMEkb <-8.62e-05# data.frame(par = names(coef(fit_all)), est = coef(fit_all)) |> # pivot_wider(names_from = par, values_from = est) |> # summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))# # get_topt# This rTPC function just finds the temperature where the function is maximized, I use the same approach brms fitsget_topt(fit_all)
[1] 9.741584
# Bootstrapping uncertainty... ? Need to refit I thinkfit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dat,start = start,lower = start*0.5,upper = start*2)# bootstrap using case resamplinglibrary(car)
Loading required package: carData
Attaching package: 'car'
The following object is masked from 'package:dplyr':
recode
The following object is masked from 'package:purrr':
some
drop_na: no rows removed
mutate: new variable 'iter' (integer) with 10,000 unique values and 0% NA
group_by_all: 5 grouping variables (r_tref, e, eh, th, iter)
ungroup: no grouping variables
mutate: new variable 'pred' (double) with 10,000,000 unique values and 0% NA
# Find the maximum for each boostrapped parameter combination ("iter")boot1_preds |>group_by(iter) |>filter(pred ==max(pred)) |>ggplot(aes(temp)) +geom_histogram()
group_by: one grouping variable (iter)
filter (grouped): removed 9,990,000 rows (>99%), 10,000 rows remaining
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Plot Sharpe Schoolfield fits
# Plot growth coefficients by year and area against mean SSTp1 <- preds |>ggplot(aes(temp, .fitted, color =factor(area, order$area))) +geom_point(data = dat, aes(temp, rate, color =factor(area, order$area)), size =1, alpha =0.8) +geom_line(linewidth =1) +geom_line(data = pred_all, aes(temp, .fitted), linewidth =1, inherit.aes =FALSE, linetype =2) +scale_color_manual(values = colors, name ="Area") +guides(color =guide_legend(nrow =1, reverse =TRUE, title.position ="top", title.hjust =0.5,override.aes =list(size =1))) +scale_x_continuous(breaks =seq(-5, 20, 2)) +labs(x ="Temperature (°C)",y ="von Bertalanffy growth coefficient (*k*)") +theme(axis.title.y = ggtext::element_markdown(),legend.position ="bottom",legend.direction ="horizontal")# Now facet and add ribbonsp1 +facet_wrap(~factor(area, levels = (arrange(area_attr, desc(lat)))$area)) +geom_ribbon(aes(ymin = lwr, ymax = upr, fill =factor(area, order$area)), alpha =0.3, color =NA) +geom_ribbon(data = pred_all, aes(temp, .fitted, ymin = lwr, ymax = upr), alpha =0.3, color =NA) +scale_fill_manual(values = colors, name ="Area") +guides(fill ="none") +geom_line(linewidth =1)
Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf
Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf
ggsave(paste0(home, "/figures/supp/sharpe_school_nls.pdf" ), width =17, height =17, units ="cm")
Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf
Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf
Can we fit a single Sharpe Schoolfield with area-specific parameters with brms?
# Again, here are the data we are fitting:ggplot(dat, aes(temp, rate, color =factor(area, order$area))) +geom_point(data = dat, aes(temp, rate, color =factor(area, order$area)), size =0.6)
# Which family to use? Well we can try with a Gaussian or a student-t first... # After that, maybe gamma or lognormal, if we predict negative K'shist(dat$rate)
# Add in fixed parametersdat$bk <-8.62e-05dat$tref <-8+273.15# We can use the nls() parameters as starting valuesavg_nls_pars <- nls_pars |>group_by(parameter) |>summarise(mean_par =mean(Estimate)) |>pivot_wider(names_from = parameter, values_from = mean_par)
group_by: one grouping variable (parameter)
summarise: now 4 rows and 2 columns, ungrouped
pivot_wider: reorganized (parameter, mean_par) into (e, eh, r_tref, th) [was 4x2, now 1x4]
avg_nls_pars
# A tibble: 1 × 4
e eh r_tref th
<dbl> <dbl> <dbl> <dbl>
1 0.611 2.75 0.309 13.0
# As a first guess, I use the nls estimate as the mean, and standard deviation something close to it. We can also bound som fo them by 0, since negative values shouldn't be possible# (better visualization including bounds further down)n=10000hist(rnorm(avg_nls_pars$e, 0.5, n=n), main ="e", xlim =c(0, 5))
hist(rnorm(avg_nls_pars$eh, 1, n=n), main ="eh", xlim =c(0, 5))
hist(rnorm(avg_nls_pars$r_tref, 0.5, n=n), main ="r_tref", xlim =c(0, 2.5))
hist(rnorm(avg_nls_pars$th, 5, n=n), main ="th", xlim =c(0, 30))
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess
# Global prior predictive check in relation to data. Doesn't loo too informative...dat |>data_grid(temp =seq_range(temp, n =51)) |>ungroup() |>mutate(bk =8.62e-05,tref =8+273.15) |>add_epred_draws(fit_prior) |>mutate(r_tref = nls_pooled_pars$r_tref, e = nls_pooled_pars$e, eh = nls_pooled_pars$eh, th = nls_pooled_pars$th) |>mutate(man_pred = (r_tref *exp(e/bk * (1/tref -1/(temp +273.15)))) / (1+exp(eh/bk * (1/(th +273.15) -1/(temp +273.15))))) |>ggplot(aes(temp, y = .epred)) +stat_lineribbon(.width =c(0.95), alpha =0.3, color ="black", fill ="black") +geom_line(aes(y = man_pred), color ="tomato3", linetype =2) +geom_point(data = dat, aes(temp, rate, color =factor(area, order$area)), size =0.6) +scale_color_manual(values = colors, name ="Area") +NULL
ungroup: no grouping variables
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
mutate (grouped): new variable 'r_tref' (double) with one unique value and 0% NA
new variable 'e' (double) with one unique value and 0% NA
new variable 'eh' (double) with one unique value and 0% NA
new variable 'th' (double) with one unique value and 0% NA
mutate (grouped): new variable 'man_pred' (double) with 51 unique values and 0% NA
# Now make sure it converges with real data but no random effectsfit_global <-brm(bf(rate ~ (rtref *exp(e/bk * (1/tref -1/(temp +273.15)))) / (1+exp(eh/bk * (1/(th +273.15) -1/(temp +273.15)))), rtref + e + eh + th ~1,nl =TRUE),data = dat,cores =4,chains =4,iter =4000, seed =9,sample_prior ="yes",prior = prior,control =list(adapt_delta =0.99, max_treedepth =12))
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
plot(fit_global)
pp_check(fit_global)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Ok, seems to work with priors and everything. They are roughly as broad as can be and still have a fit. Now fit the full model, with random area effects, more iterations and chains!
# Still good enough, might change some priors. Now look at random area effects!# Gaussian modelfitb <-brm(bf(rate ~ (rtref *exp(e/bk * (1/tref -1/(temp +273.15)))) / (1+exp(eh/bk * (1/(th +273.15) -1/(temp +273.15)))), rtref + e + eh + th ~1+ (1|area),nl =TRUE),data = dat,cores =4,chains =4,iter =5000,sample_prior ="yes",save_pars =save_pars(all =TRUE),seed =9,prior = prior,control =list(adapt_delta =0.99, max_treedepth =12))
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
print(fitbs, digits =4)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))
rtref ~ 1 + (1 | area)
e ~ 1 + (1 | area)
eh ~ 1 + (1 | area)
th ~ 1 + (1 | area)
Data: dat (Number of observations: 306)
Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup draws = 10000
Group-Level Effects:
~area (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(rtref_Intercept) 0.0322 0.0232 0.0020 0.0899 1.0007 1732
sd(e_Intercept) 0.1879 0.1283 0.0109 0.4904 1.0008 2638
sd(eh_Intercept) 0.4462 0.4463 0.0150 1.5712 1.0006 2023
sd(th_Intercept) 1.1362 1.0236 0.0454 3.6889 1.0002 3409
Tail_ESS
sd(rtref_Intercept) 2690
sd(e_Intercept) 3115
sd(eh_Intercept) 4198
sd(th_Intercept) 4052
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept 0.3152 0.1168 0.1962 0.6317 1.0012 1404 2871
e_Intercept 0.5867 0.3309 0.0879 1.3781 1.0010 1627 2523
eh_Intercept 1.7546 0.6191 0.7373 3.3783 1.0025 2292 1961
th_Intercept 12.2022 4.1609 4.7838 20.7330 1.0011 1529 2959
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.0397 0.0021 0.0357 0.0438 1.0009 9550 7205
nu 24.7640 13.4838 7.9340 59.1615 1.0002 9928 8129
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Compare fitted models based on ELPD. Preferred model in the first row!
fitb_loo <-loo(fitb, moment_match =TRUE)
Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
print(fitbs_flat, digits =4)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))
rtref ~ 1 + (1 | area)
e ~ 1 + (1 | area)
eh ~ 1 + (1 | area)
th ~ 1 + (1 | area)
Data: dat (Number of observations: 306)
Draws: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~area (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(rtref_Intercept) 0.0395 0.0313 0.0020 0.1182 1.0005 2315
sd(e_Intercept) 0.1880 0.1311 0.0096 0.5056 1.0006 3308
sd(eh_Intercept) 0.3715 0.3886 0.0141 1.3542 1.0031 2393
sd(th_Intercept) 1.0269 0.9518 0.0315 3.4494 1.0012 3386
Tail_ESS
sd(rtref_Intercept) 4256
sd(e_Intercept) 4662
sd(eh_Intercept) 4541
sd(th_Intercept) 4968
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept 0.3849 0.1621 0.1994 0.8170 1.0035 1642 2433
e_Intercept 0.6985 0.4325 0.0862 1.7128 1.0031 1730 4029
eh_Intercept 1.5619 0.6696 0.3278 3.2472 1.0010 3031 2471
th_Intercept 10.4526 4.4138 3.7388 19.1941 1.0039 1622 3638
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.0397 0.0020 0.0358 0.0437 1.0001 11322 8949
nu 24.7114 13.5015 8.0481 59.5055 1.0001 12506 10073
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
print(fitbs, digits =4)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))
rtref ~ 1 + (1 | area)
e ~ 1 + (1 | area)
eh ~ 1 + (1 | area)
th ~ 1 + (1 | area)
Data: dat (Number of observations: 306)
Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup draws = 10000
Group-Level Effects:
~area (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(rtref_Intercept) 0.0322 0.0232 0.0020 0.0899 1.0007 1732
sd(e_Intercept) 0.1879 0.1283 0.0109 0.4904 1.0008 2638
sd(eh_Intercept) 0.4462 0.4463 0.0150 1.5712 1.0006 2023
sd(th_Intercept) 1.1362 1.0236 0.0454 3.6889 1.0002 3409
Tail_ESS
sd(rtref_Intercept) 2690
sd(e_Intercept) 3115
sd(eh_Intercept) 4198
sd(th_Intercept) 4052
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept 0.3152 0.1168 0.1962 0.6317 1.0012 1404 2871
e_Intercept 0.5867 0.3309 0.0879 1.3781 1.0010 1627 2523
eh_Intercept 1.7546 0.6191 0.7373 3.3783 1.0025 2292 1961
th_Intercept 12.2022 4.1609 4.7838 20.7330 1.0011 1529 2959
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.0397 0.0021 0.0357 0.0438 1.0009 9550 7205
nu 24.7640 13.4838 7.9340 59.1615 1.0002 9928 8129
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
The student model is preferred. Inspect it!
# Check fitsummary(fitbs)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))
rtref ~ 1 + (1 | area)
e ~ 1 + (1 | area)
eh ~ 1 + (1 | area)
th ~ 1 + (1 | area)
Data: dat (Number of observations: 306)
Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup draws = 10000
Group-Level Effects:
~area (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(rtref_Intercept) 0.03 0.02 0.00 0.09 1.00 1732 2690
sd(e_Intercept) 0.19 0.13 0.01 0.49 1.00 2638 3115
sd(eh_Intercept) 0.45 0.45 0.02 1.57 1.00 2023 4198
sd(th_Intercept) 1.14 1.02 0.05 3.69 1.00 3409 4052
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept 0.32 0.12 0.20 0.63 1.00 1404 2871
e_Intercept 0.59 0.33 0.09 1.38 1.00 1627 2523
eh_Intercept 1.75 0.62 0.74 3.38 1.00 2292 1961
th_Intercept 12.20 4.16 4.78 20.73 1.00 1529 2959
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.04 0.00 0.04 0.04 1.00 9550 7205
nu 24.76 13.48 7.93 59.16 1.00 9928 8129
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
ts <- dat |>data_grid(temp =seq_range(temp, n =100)) |>mutate(bk =8.62e-05,tref =8+273.15) |>add_epred_draws(fitbs, re_formula =NA, ndraws =10000)
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ggplot(ts) +geom_line(aes(temp, y = .epred, group = .draw), alpha = .1)
# Compute quantiles across spaghettiest_opt <- ts |>group_by(.draw) |>filter(.epred ==max(.epred)) |>ungroup() |>filter(!temp ==min(temp)) |># remove values where there was no optimumfilter(!temp ==max(temp))
# To add the histogram inside the main plot, I can convert it to "numbers-at-temperature" and plot as a bin. This way I know the maximum X. I will plot these as separatly. In ggplot, you cannot technically do this, but you can scale the axes to have similar values, the change labels.t_opt2 <- t_opt |>group_by(temp) |>summarise(n =n())
group_by: one grouping variable (temp)
summarise: now 98 rows and 2 columns, ungrouped
ggplot(t_opt2, aes(temp, n)) +geom_bar(stat ="identity", color =NA, fill ="grey30", width =0.14)
Warning: `position_stack()` requires non-overlapping x intervals
scaleFactor <-max(t_opt2$n) /max(dat$rate)# Make the main plot (conditional effect of temperature, with and without random effects)# Predictions without random effectsglob <- dat |>data_grid(temp =seq_range(temp, n =100)) |>mutate(bk =8.62e-05,tref =8+273.15) |>add_epred_draws(fitbs, re_formula =NA) |>ungroup()
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
group_by: one grouping variable (area)
ungroup: no grouping variables
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
Warning: `position_stack()` requires non-overlapping x intervals
ggsave(paste0(home, "/figures/sharpe_school_bayes.pdf" ), width =17, height =17, units ="cm")
Warning: `position_stack()` requires non-overlapping x intervals
Area-specific T_opts as a ridgeplot
# Calculate T_opt by areatsa <- dat |>data_grid(area =unique(dat$area),temp =seq_range(temp, n =100)) |>mutate(bk =8.62e-05,tref =8+273.15) |>add_epred_draws(fitbs, re_formula =NULL, ndraws =10000)
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
t_opt_a <- tsa |>group_by(.draw, area) |>filter(.epred ==max(.epred)) |>ungroup() |>filter(!temp ==min(temp)) |># remove values where there was no optimumfilter(!temp ==max(temp))
mutate (grouped): changed 400,000 values (100%) of 'parameter' (0 new NA)
Adding missing grouping variables: `Intercept`
left_join: added one column (.pop_value)
> rows only in x 0
> rows only in y ( 0)
> matched rows 400,000
> =========
> rows total 400,000
mutate (grouped): changed 400,000 values (100%) of '.value' (0 new NA)
summarise: now 40 rows and 6 columns, 2 group variables remaining (area,
Intercept)
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA
pars <-bind_rows(brms_pars, nls_pars)# Need to trim errorbars...ggplot(pars, aes(area, Estimate, color = source)) +geom_point() +geom_errorbar(aes(ymin = lwr, ymax = upr), width =0) +facet_wrap(~parameter, scales ="free")
# And without errorbars...pars |>ggplot(aes(area, Estimate, color = source)) +geom_point(position =position_dodge(width =0.5)) +facet_wrap(~parameter, scales ="free")
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
rename: renamed one variable (parameter)
mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)
summarise: now 4 rows and 4 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA
new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 289 rows (94%), 17 rows remaining
Compiling Stan program...
recompiling to avoid crashing R session
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
rename: renamed one variable (parameter)
mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)
summarise: now 4 rows and 4 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA
new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 259 rows (85%), 47 rows remaining
Compiling Stan program...
recompiling to avoid crashing R session
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
rename: renamed one variable (parameter)
mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)
summarise: now 4 rows and 4 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA
new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 269 rows (88%), 37 rows remaining
Compiling Stan program...
Start sampling
rename: renamed one variable (parameter)
mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)
summarise: now 4 rows and 4 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA
new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 277 rows (91%), 29 rows remaining
Compiling Stan program...
Start sampling
rename: renamed one variable (parameter)
mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)
summarise: now 4 rows and 4 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA
new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 246 rows (80%), 60 rows remaining
Compiling Stan program...
Start sampling
rename: renamed one variable (parameter)
mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)
summarise: now 4 rows and 4 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA
new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 288 rows (94%), 18 rows remaining
Compiling Stan program...
Start sampling
rename: renamed one variable (parameter)
mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)
summarise: now 4 rows and 4 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA
new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 292 rows (95%), 14 rows remaining
Compiling Stan program...
Start sampling
rename: renamed one variable (parameter)
mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)
summarise: now 4 rows and 4 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA
new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 268 rows (88%), 38 rows remaining
Compiling Stan program...
Start sampling
rename: renamed one variable (parameter)
mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)
summarise: now 4 rows and 4 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA
new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 282 rows (92%), 24 rows remaining
Compiling Stan program...
Start sampling
rename: renamed one variable (parameter)
mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)
summarise: now 4 rows and 4 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA
new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
mutate: new variable 'area' (character) with one unique value and 0% NA
# Save aic as data framesarea_spec <-bind_rows(dlist)area_spec_preds <-bind_rows(dpred)# Plot prediction comparisonggplot(area_pred_brms, aes(temp, y = .epred, color =factor(area, order$area),fill =factor(area, order$area))) +stat_lineribbon(.width =c(0.95), alpha =0.1, color =NA) +stat_lineribbon(.width =c(0), alpha =0.6) +# add non-hierarchical model predictionsstat_lineribbon(data = area_spec_preds, .width =0, alpha =0.6,aes(temp, y = .epred, color =factor(area, order$area), linetype ="non-hierarchial")) +geom_point(data = dat, aes(temp, rate, color =factor(area, order$area)), size =0.4, alpha =0.4) +scale_color_manual(values = colors, name ="Area") +scale_fill_manual(values = colors, name ="Area") +scale_linetype_manual(values =2) +facet_wrap(~factor(area, rev(order$area))) +guides(fill ="none",color =guide_legend(nrow =1, reverse =TRUE, title.position ="top", title.hjust =0.5,override.aes =list(size =1))) +theme(axis.title.y = ggtext::element_markdown(),legend.position ="bottom",legend.direction ="horizontal") +labs(x ="Temperature (°C)",y ="von Bertalanffy growth coefficient (*k*)",linetype ="")
# Now zoom in on RA and add the nls fitggplot(area_pred_brms |>filter(area =="RA"), aes(temp, y = .epred, color =factor(area, order$area),fill =factor(area, order$area))) +stat_lineribbon(.width =c(0.95), alpha =0.1, color =NA) +stat_lineribbon(.width =c(0), alpha =0.6) +# add non-hierarchical model predictionsstat_lineribbon(data = area_spec_preds |>filter(area =="RA"), .width =0, alpha =0.6,aes(temp, y = .epred, color =factor(area, order$area), linetype ="non-hierarchial")) +# add nlsgeom_line(data = preds |>filter(area =="RA"), aes(temp, .fitted, linetype ="nls multstart"), color ="gray10") +geom_point(data = dat |>filter(area =="RA"), aes(temp, rate, color =factor(area, order$area)), size =0.4, alpha =0.4) +scale_color_manual(values = colors, name ="Area") +scale_fill_manual(values = colors, name ="Area") +scale_linetype_manual(values =c(2, 3)) +facet_wrap(~factor(area, rev(order$area))) +guides(fill ="none",color =guide_legend(nrow =1, reverse =TRUE, title.position ="top", title.hjust =0.5,override.aes =list(size =1))) +theme(axis.title.y = ggtext::element_markdown(),legend.position ="bottom",legend.direction ="horizontal") +labs(x ="Temperature (°C)",y ="von Bertalanffy growth coefficient (*k*)",linetype ="")
Temperature in growing season instead of all year average
Fit Sharpe-Schoolfield model to K
By area
model <-'sharpeschoolhigh_1981'# Get starting values on full dataset for Sharpe-Schoolfield modeldat2 <- VBG |>select(k_median, temp_gs, area) |>rename(rate = k_median,temp = temp_gs) # growing season! not annual average!lower <-get_lower_lims(dat2$temp, dat2$rate, model_name = model)upper <-get_upper_lims(dat2$temp, dat2$rate, model_name = model)start <-get_start_vals(dat2$temp, dat2$rate, model_name = model)# Sharpe-Schoolfield model fit to data for each areapreds2 <-NULLpred2 <-NULLfor(a inunique(dat2$area)) {# Get data dd <- dat2 |>filter(area == a)# Fit model fit2 <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dd,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new data new_data2 <-data.frame(temp =seq(min(dd$temp), max(dd$temp), length.out =100)) pred2 <-augment(fit2, newdata = new_data2) |>mutate(area = a)# Add to general data frame preds2 <-data.frame(rbind(preds2, pred2))}
All areas pooled
# One sharpe schoolfield fitted to all datafit_all2 <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dat2,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new datanew_data_all2 <-data.frame(temp =seq(min(dat2$temp), max(dat2$temp), length.out =100))pred_all2 <-augment(fit_all2, newdata = new_data_all2) |>mutate(area ="all")# Add t_optkb <-8.62e-05data.frame(par =names(coef(fit_all2)), est =coef(fit_all2)) |>pivot_wider(names_from = par, values_from = est) |>summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))
Plot Sharpe Schoolfield fits
# Plot growth coefficients by year and area against mean SSTpreds2 |>ggplot(aes(temp, .fitted, color =factor(area, order$area))) +geom_point(data = dat2, aes(temp, rate, color =factor(area, order$area)), size =0.6) +geom_line(linewidth =1) +geom_line(data = pred_all2, aes(temp, .fitted), linewidth =1, inherit.aes =FALSE, linetype =2) +scale_color_manual(values = colors, name ="Area") +guides(color =guide_legend(override.aes =list(size =1))) +scale_x_continuous(breaks =seq(-5, 20, 1)) +labs(x ="Temperature (°C)",y ="von Bertalanffy growth coefficient (*k*)") +theme(axis.title.y = ggtext::element_markdown())
Linf instead of k
By area
model <-'sharpeschoolhigh_1981'# Get starting values on full dataset for Sharpe-Schoolfield modeldat3 <- VBG |>select(linf_median, temp, area) |>rename(rate = linf_median)lower <-get_lower_lims(dat3$temp, dat3$rate, model_name = model)upper <-get_upper_lims(dat3$temp, dat3$rate, model_name = model)start <-get_start_vals(dat3$temp, dat3$rate, model_name = model)# Sharpe-Schoolfield model fit to data for each areapreds3 <-NULLpred3 <-NULLfor(a inunique(dat3$area)) {# Get data dd <- dat3 |>filter(area == a)# Fit model fit3 <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dd,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new data new_data3 <-data.frame(temp =seq(min(dd$temp), max(dd$temp), length.out =100)) pred3 <-augment(fit3, newdata = new_data3) |>mutate(area = a)# Add to general data frame preds3 <-data.frame(rbind(preds3, pred3))}
All areas pooled
# One sharpe schoolfield fitted to all datafit_all3 <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dat3,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new datanew_data_all3 <-data.frame(temp =seq(min(dat3$temp), max(dat3$temp), length.out =100))pred_all3 <-augment(fit_all3, newdata = new_data_all3) |>mutate(area ="all")
Plot Sharpe Schoolfield fits
# Plot growth coefficients by year and area against mean SSTpreds3 |>ggplot(aes(temp, .fitted, color =factor(area, order$area))) +geom_point(data = dat3, aes(temp, rate, color =factor(area, order$area)), size =0.6) +geom_line(linewidth =1) +geom_line(data = pred_all3, aes(temp, .fitted), linewidth =1, inherit.aes =FALSE, linetype =2) +scale_color_manual(values = colors, name ="Area") +guides(color =guide_legend(override.aes =list(size =1))) +scale_x_continuous(breaks =seq(-5, 20, 1)) +labs(x ="Temperature (°C)",y =expression(paste(italic(L[infinity]), " [cm]"))) +theme(axis.title.y = ggtext::element_markdown())